# packages
require(arm)
require(data.table)
require(effects)
require(foreach)
require(here)
require(ggimage)
require(ggplot2)
require(ggpubr)
require(grid)
require(gtable)
require(MASS)
require(multcomp)
require(optimx)
require(performance)
require(PerformanceAnalytics)
require(png)
require(RColorBrewer)
require(rphylopic)
require(viridis)
# constants
round_ = 3 # number of decimal places to round model coefficients
nsim = 5000 # number of simulations to extract estimates and 95%CrI
ax_lines = "grey60" # defines color of the axis lines
#colors <- c("#999999", "#E69F00", "#56B4E9") #viridis(3)
set.seed(42)
width_ = .7 # spacing between error bars
col_ = c(brewer.pal(n =12, name = "Paired"), 'grey30','grey80')
# functions
# to add images to panels
annotation_custom2 <- function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data) {
layer(data = data, stat = StatIdentity, position = PositionIdentity,
geom = ggplot2:::GeomCustomAnn,
inherit.aes = TRUE, params = list(grob = grob,
xmin = xmin, xmax = xmax,
ymin = ymin, ymax = ymax))
}
# to remove ggplot components
gtable_filter_remove <- function (x, name, trim = TRUE){
matches <- !(x$layout$name %in% name)
x$layout <- x$layout[matches, , drop = FALSE]
x$grobs <- x$grobs[matches]
if (trim)
x <- gtable_trim(x)
x
}
# customized ggplot theme
theme_MB = theme(
title = element_text(size=8, colour="grey30"),
axis.line = element_blank(),
#axis.line = element_line(colour="grey70", size=0.25),
axis.title = element_text(size=7, colour="grey30"),
axis.title.y = element_text(vjust=3.5),
axis.title.x = element_text(vjust=1),
axis.text = element_text(size=6),#, vjust = 0.5, hjust=1),# margin=units(0.5,"mm")),
axis.ticks.length=unit(0.5,"mm"),
axis.ticks = element_line(colour = "grey70", size = 0.1),
#axis.ticks.margin,
strip.text.x = element_text(size = 6, color="grey30", margin=margin(1,1,1,1,"mm")), #grey50
strip.text.y = element_text(size = 6, color="grey30", margin=margin(1,1,1,1,"mm")), #grey50
strip.background = element_rect(fill="grey99",colour="grey70", size=0.25),
#strip.background = element_blank(),
#strip.text = element_blank(),
panel.spacing = unit(0, "mm"),
panel.background=element_blank(),
panel.border = element_rect(colour="grey70", size=0.1, fill = NA), #panel.border=element_blank(),
panel.grid = element_blank(),
legend.text=element_text(size=6),
legend.title=element_text(size=6),
legend.key = element_rect(colour = NA, fill = NA),
legend.key.height= unit(0.5,"line"),
legend.key.width = unit(0.25, "cm"),
legend.margin = margin(0,0,0,0, unit="cm"),
legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
legend.background = element_blank()
)
# for estimates
est_out =function(model = m, label = "", nsim = 5000){
bsim = sim(model, n.sim=5000)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975))
data.table(predictor=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,], model = paste(label, "N =", nobs(model)))
}
# change color
change_col = function(replace_black, theimg) {
r_b = col2rgb(replace_black) / 255
#theimg[theimg == 1] <- 2
for (i in 1:3) {
theimg[,,i][theimg[,,i] == 0] <- r_b[i]
}
return(theimg)
}
# for Supplementary Table output based on sim
m_out = function(model = m, type = "mixed",
name = "define", dep = "define", fam = 'Gaussian',
round_ = 3, nsim = 5000, aic = FALSE, save_sim = 'Data/model_sim/', back_tran = FALSE, perc_ = 1){
# perc_ 1 = proportion or 100%
bsim = sim(model, n.sim=nsim)
if(save_sim!=FALSE){save(bsim, file = paste0(save_sim, name,'.RData'))}
if(type != "mixed"){
v = apply(bsim@coef, 2, quantile, prob=c(0.5))
ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975))
if(back_tran == TRUE & fam == "binomial"){
v = perc_*plogis(v)
ci = perc_*plogis(ci)
}
if(back_tran == TRUE & fam == "binomial_logExp"){
v = perc_*(1-plogis(v))
ci = perc_*(1-plogis(ci))
ci = rbind(ci[2,],ci[1,])
}
if(back_tran == TRUE & fam == "Poisson"){
v = exp(v)
ci = exp(ci)
}
oi=data.frame(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
rownames(oi) = NULL
oi$estimate_r=round(oi$estimate,round_)
oi$lwr_r=round(oi$lwr,round_)
oi$upr_r=round(oi$upr,round_)
if(perc_ == 100){
oi$estimate_r = paste0(oi$estimate_r,"%")
oi$lwr_r = paste0(oi$lwr_r,"%")
oi$upr_r = paste0(oi$upr_r,"%")
}
x=data.table(oi[c('type',"effect", "estimate_r","lwr_r",'upr_r')])
}else{
v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975))
if(back_tran == TRUE & fam == "binomial"){
v = perc_*plogis(v)
ci = perc_*plogis(ci)
}
if(back_tran == TRUE & fam == "binomial_logExp"){
v = perc_*(1-plogis(v))
ci = perc_*(1-plogis(ci))
ci = rbind(ci[2,],ci[1,])
}
if(back_tran == TRUE & fam == "Poisson"){
v = exp(v)
ci = exp(ci)
}
oi=data.table(type='fixed',effect=rownames(coef(summary(model))),estimate=v, lwr=ci[1,], upr=ci[2,])
rownames(oi) = NULL
oi[,estimate_r := round(estimate,round_)]
oi[,lwr_r := round(lwr,round_)]
oi[,upr_r :=round(upr,round_)]
if(perc_ == 100){
oi[,estimate_r := paste0(estimate_r,"%")]
oi[,lwr_r := paste0(lwr_r,"%")]
oi[,upr_r := paste0(upr_r,"%")]
}
oii=oi[,c('type',"effect", "estimate_r","lwr_r",'upr_r')]
l=data.frame(summary(model)$varcor)
l = l[is.na(l$var2),]
l$var1 = ifelse(is.na(l$var1),"",l$var1)
l$pred = paste(l$grp,l$var1)
q050={}
q025={}
q975={}
pred={}
# variance of random effects
for (ran in names(bsim@ranef)) {
ran_type = l$var1[l$grp == ran]
for(i in ran_type){
q050=c(q050,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.5)))
q025=c(q025,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.025)))
q975=c(q975,quantile(apply(bsim@ranef[[ran]][,,ran_type], 1, var), prob=c(0.975)))
pred= c(pred,paste(ran, i))
}
}
# residual variance
q050=c(q050,quantile(bsim@sigma^2, prob=c(0.5)))
q025=c(q025,quantile(bsim@sigma^2, prob=c(0.025)))
q975=c(q975,quantile(bsim@sigma^2, prob=c(0.975)))
pred= c(pred,'Residual')
ri=data.table(type='random %',effect=pred, estimate_r=round(100*q050/sum(q050)), lwr_r=round(100*q025/sum(q025)), upr_r=round(100*q975/sum(q975)))
ri[lwr_r>upr_r, lwr_rt := upr_r]
ri[lwr_r>upr_r, upr_rt := lwr_r]
ri[!is.na(lwr_rt), lwr_r := lwr_rt]
ri[!is.na(upr_rt), upr_r := upr_rt]
ri$lwr_rt = ri$upr_rt = NULL
ri[,estimate_r := paste0(estimate_r,'%')]
ri[,lwr_r := paste0(lwr_r,'%')]
ri[,upr_r := paste0(upr_r,'%')]
x = data.table(rbind(oii,ri))
}
x[1, model := name]
x[1, response := dep]
x[1, error_structure := fam]
N = length(resid(model))
x[1, N := N ]
x=x[ , c('model', 'response', 'error_structure', 'N', 'type',"effect", "estimate_r","lwr_r",'upr_r')]
if (aic == TRUE){
x[1, AIC := AIC(update(model,REML = FALSE))]
}
if (aic == "AICc"){
aicc = AICc(model)
x[1, AICc := aicc]
}
if(type == "mixed" & nrow(x[type=='random %' & estimate_r =='0%'])==0){
R2_mar = as.numeric(r2_nakagawa(model)$R2_marginal)
R2_con = as.numeric(r2_nakagawa(model)$R2_conditional)
x[1, R2_mar := R2_mar]
x[1, R2_con := R2_con]
}
x[is.na(x)] = ""
return(x)
}
# model assumption function
m_ass = function(name = 'define', mo = m0, dat = d, fixed = NULL, categ = NULL, trans = "none", spatial = TRUE, temporal = TRUE, PNG = TRUE, outdir = 'outdir', n_col=8, width_ = 10, height_ = 5.5){
l=data.frame(summary(mo)$varcor)
l = l[is.na(l$var2),]
nt = if(temporal==TRUE){1}else{0}
ns = if(spatial==TRUE){7}else{0}
n = 3+nrow(l)-1+length(fixed)+length(categ) + nt + ns
if(PNG == TRUE){
png(paste(outdir,name, ".png", sep=""), width=width_,height=height_,units="in",res=600) # width = 6
par(mfrow=c(4, n_col),tcl = -0.08, cex = 0.5, cex.main = 0.9,#ceiling(n/n_col),n_col)
oma = c(1,1,2,1),
mar = c(2, 2, 2, 1), mgp=c(1,0,0)
)
}else{
dev.new(width=width_,height=height_)
par(mfrow=c(4,n_col), tcl = -0.08, cex = 0.5, cex.main = 0.9,#ceiling(n/n_col),n_col)
oma = c(1,1,2,1),
mar = c(2, 2, 2, 1), mgp=c(1,0,0)
)
}
scatter.smooth(fitted(mo),resid(mo),col='grey');abline(h=0, lty=2, col ='red')
scatter.smooth(fitted(mo),sqrt(abs(resid(mo))), col='grey')
qqnorm(resid(mo), main=list("Normal Q-Q Plot: residuals"),col='grey');qqline(resid(mo), col = 'red')
#unique(l$grp[l$grp!="Residual"])
for(i in unique(l$grp[l$grp!="Residual"])){
#i = "mean_year"
ll=ranef(mo)[names(ranef(mo))==i][[1]]
if(ncol(ll)==1){
qqnorm(ll[,1], main = paste(i,names(ll)[1]),col='grey',);qqline(ll[,1], col ='red')
}else{
qqnorm(ll[,1], main = paste(i,names(ll)[1]),col='grey');qqline(ll[,1], col ='red')
qqnorm(ll[,2], main = paste(i,names(ll)[2]),col='grey');qqline(ll[,2], col ='red')
}
}
# variables
scatter={}
for (i in rownames(summary(mo)$coef)) {
#i = "lat_abs"
j=sub("\\).*", "", sub(".*\\(", "",i))
scatter[length(scatter)+1]=j
}
x = data.frame(scatter=unique(scatter)[2:length(unique(scatter))],
log_ = grepl("log",rownames(summary(mo)$coef)[2:length(unique(scatter))]), stringsAsFactors = FALSE)
for (i in 1:length(fixed)){
jj =fixed[i]
variable=dat[, ..jj][[1]]
if(trans[i]=='log'){
scatter.smooth(resid(mo)~log(variable),xlab=paste('log(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}else if(trans[i]=='abs'){
scatter.smooth(resid(mo)~abs(variable),xlab=paste('abs(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}else if(trans[i]=='sin'){scatter.smooth(resid(mo)~sin(variable),xlab=paste('sin(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}else if(trans[i]=='cos'){scatter.smooth(resid(mo)~cos(variable),xlab=paste('cos(',jj,')',sep=''), col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}else{
scatter.smooth(resid(mo)~variable,xlab=jj,col = 'grey');abline(h=0, lwd=1, lty = 2, col ='red')
}
}
if(length(categ)>0){
for(i in categ){
variable=dat[, ..i][[1]]
boxplot(resid(mo)~variable, medcol='grey', whiskcol='grey', staplecol='grey', boxcol='grey', outcol='grey', xlab = i);abline(h=0, lty=3, lwd=1, col = 'red')
}
}
if(temporal == TRUE){
acf(resid(mo), type="p", main=list("Temporal autocorrelation:\npartial series residual"))
}
if(spatial == TRUE){
spdata=data.frame(resid=resid(mo), x=dat$Lon, y=dat$Lat)
spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
#cex_=c(1,2,3,3.5,4)
cex_=c(1,1.5,2,2.5,3)
spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
plot(spdata$x, spdata$y,col=spdata$col, cex=as.numeric(spdata$cex), pch= 16, main=list('Spatial distribution of residuals', cex=0.8), xlab = 'Longitude', ylab = 'Latitude')
legend("topright", pch=16, legend=c('>0','<0'), ,col=c(rgb(83,95,124,100, maxColorValue = 255),rgb(253,184,19,100, maxColorValue = 255)))
plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('residual >=0'), xlab = 'Longitude', ylab = 'Latitude')
# EU
dat$res = resid(mo)
spdata=data.frame(resid = dat$res[dat$Country!='Australia'], x=dat$Lon[dat$Country!='Australia'], y=dat$Lat[dat$Country!='Australia'])
spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
#cex_=c(1,2,3,3.5,4)
cex_=c(1,1.5,2,2.5,3)
spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('EU - residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('EU residuals >=0)'), xlab = 'Longitude', ylab = 'Latitude')
# Australia
spdata=data.frame(resid = dat$res[dat$Country=='Australia'], x=dat$Lon[dat$Country=='Australia'], y=dat$Lat[dat$Country=='Australia'])
spdata$col=ifelse(spdata$resid<0,rgb(83,95,124,100, maxColorValue = 255),ifelse(spdata$resid>0,rgb(253,184,19,100, maxColorValue = 255), 'red'))
#cex_=c(1,2,3,3.5,4)
cex_=c(1,1.5,2,2.5,3)
spdata$cex=as.character(cut(abs(spdata$resid), 5, labels=cex_))
plot(spdata$x[spdata$resid<0], spdata$y[spdata$resid<0],col=spdata$col[spdata$resid<0], cex=as.numeric(spdata$cex[spdata$resid<0]), pch= 16, main=list('Australia residuals <0'), xlab = 'Longitude', ylab = 'Latitude')
plot(spdata$x[spdata$resid>=0], spdata$y[spdata$resid>=0],col=spdata$col[spdata$resid>=0], cex=as.numeric(spdata$cex[spdata$resid>=0]), pch= 16, main=list('Australia residuals >=0'), xlab = 'Longitude', ylab = 'Latitude')
}
mtext(stringr::str_wrap(paste(paste0(name," model: "), slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''), width = ceiling(nchar(paste(slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''))/2)+10), side = 3, line = 0, cex=0.5,outer = TRUE, col = 'darkblue') #ceiling(nchar(paste(slot(mo,"call")[1],'(',slot(mo,"call")[2],sep=''))/2)
if(PNG==TRUE){dev.off()}
}
# data
o = fread(here::here('Data/phylopic.txt'))
setnames(o, old = c('Name', 'Code'), new = c('genus2', 'uid'))
t = fread(here::here('Data/taxonomy.txt'))
g = fread(here::here('Data/google_mobility.txt')) #fwrite(d, here::here('Data/data.txt'), sep ='\t')
g[country_region == "Czechia", country_region := 'Czech Republic']
g[, Year := as.integer(substring(date, nchar(date)-3, nchar(date)))]
g[nchar(date)==9, date:=paste0('0',date)]
g[, date_ :=as.POSIXct(date, format = '%d.%m.%y')]
g[, Day :=yday(date_)]
g[country_region!='Ausstralia', Day := Day-92 +1] # 1 April = start of breeding season (1st day) = 92 day of the year
g[country_region=='Ausstralia', Day := Day-228 +1] # 15 Augusst = start of breeding season (1st day) = 228 day of the year
setnames(g, old = 'country_region', new ='Country')
d = fread(here::here('Data/data_corrected.txt')) #fwrite(d, here::here('Data/data.txt'), sep ='\t')
# adjust correct assignment of season (Year) for Australia
d[Country == 'Australia' & Year == 2020 & Covid == 0, Year:=2019]
d[Country == 'Australia' & Year == 2021 & Day>139, Year:=2020]
d[Country == 'Australia' & Year == 2022 & Day>139, Year:=2021]
d = d[order(Year, IDLocality, Day, Hour)]
d[Country%in%'Czech_Republic', Country:='Czech Republic']
d[, genus := sub("_.*", "", Species)]
d[, sp_day_year := paste(Year, Species, Day, sep="_")]
d[, sp_loc := paste(Species, IDLocality, sep="_")]
d[, sp_country := paste(Species, Country, sep="\n")]
d[, rad:=(2*pi*Hour) / 24]
d[, Day_:= Day]
#d[, FID_z := scale(FID), by = Species]
#d[, SD_z := scale(SD), by = Species]
d[, FID_ln := log(FID)]
d[, SD_ln := log(SD)]
d[, body_ln := log(BodyMass)]
d[, flock_ln := log(FlockSize)]
#d[Country == 'Australia', Day_:= abs(Day - 189)]
d1 = d[Covid == 1, .N, by = Species]
d2 = d[Covid == 0, .N, by = Species]
setnames(d1, old = 'N', new ='N_during')
setnames(d2, old = 'N', new ='N_before')
dd = merge(d1,d2) # species with data before and during
da = merge(d1,d2, all = TRUE)
d1p = d[Country!='Poland' & Covid == 1, .N, by = Species]
d2p = d[Country!='Poland' & Covid == 0, .N, by = Species]
setnames(d1p, old = 'N', new ='N_during')
setnames(d2p, old = 'N', new ='N_before')
ddp = merge(d1p,d2p) # species with data before and during, but without Poland data
p1 = d[Covid == 1, .N, by = .(IDLocality, Species)]
p2 = d[Covid == 0, .N, by = .(IDLocality, Species)]
setnames(p1, old = 'N', new ='N_during')
setnames(p2, old = 'N', new ='N_before')
pp = merge(p1,p2) # species-localities with data before and during
pa = merge(p1,p2, all = TRUE)
p1p = d[Country!='Poland' & Covid == 1, .N, by = .(IDLocality, Species)]
p2p = d[Country!='Poland' & Covid == 0, .N, by = .(IDLocality, Species)]
setnames(p1p, old = 'N', new ='N_during')
setnames(p2p, old = 'N', new ='N_before')
ppp = merge(p1p,p2p) # species-localities with data before and during,but without Poland data
p1p4 = d[Year!=2014 & Country!='Poland' & Covid == 1, .N, by = .(IDLocality, Species)]
p2p4 = d[Year!=2014 &Country!='Poland' & Covid == 0, .N, by = .(IDLocality, Species)]
setnames(p1p4, old = 'N', new ='N_during')
setnames(p2p4, old = 'N', new ='N_before')
ppp4 = merge(p1p4,p2p4) # species-localities with data before and during, but without Poland & 2014 data
s = d[Covid == 1]
s[, Nsp := .N, by ='Species']
s[, sp := gsub('[_]', ' ', Species)]
# add google mobility
s = merge(s, g[,.(Country, Year, date_,Day, parks_percent_change_from_baseline)], all.x = TRUE)
s[, year_ := as.character(Year)]
s[, weekday := weekdays(date_)]
ss = s[!is.na(parks_percent_change_from_baseline)]
ss[, country_year := paste(Country, Year)] #table(paste(s$Country, s$Year))
Google Mobility Index uses February baseline values for each city and day of the week and reports % changes for each day of the week. In other words, a 10% increase on Monday may mean differnt human mobility than 10% increase on Sunday and this may differ between cities. We can control for this in the models but it is worth keeping in mind. ## Distributions’
ggplot(g, aes(x = parks_percent_change_from_baseline, fill = factor(Year))) +
geom_histogram(position = "dodge") +
#scale_y_continuous(trans = 'log')+
scale_fill_manual(values = c('orange', 'skyblue', 'black')) +
geom_vline(xintercept = 0, lty = 3, col = 'red')+
facet_wrap(~Country, nrow = 5)
ggplot(g, aes(x = Country, y = parks_percent_change_from_baseline, fill = factor(Year))) +
geom_boxplot()+
#scale_y_continuous(trans = 'log')+
scale_fill_manual(values = c('orange', 'skyblue', 'black'))
ggplot(g, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
geom_line()+
facet_wrap(~Country, nrow = 5) +
#scale_y_continuous(trans = 'log')+
scale_color_manual(values = c('orange', 'skyblue', 'black'))
ggplot(g, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
stat_smooth()+
facet_wrap(~Country, nrow = 5) +
labs(subtitle = 'Loess' )+
#scale_y_continuous(trans = 'log')+
scale_color_manual(values = c('orange', 'skyblue', 'black'))
ggplot(g, aes(x = Day, y = parks_percent_change_from_baseline, col = factor(Year))) +
stat_smooth(method = 'lm')+
facet_wrap(~Country, nrow = 5) +
labs(subtitle = 'linear fit' )+
#scale_y_continuous(trans = 'log')+
scale_color_manual(values = c('orange', 'skyblue', 'black'))
!! Although histograms show slight shift (perhaps decline) in the distribution of human mobility in 2020 (21) in compariison to post covid 2022 and boxplots perhaps confirm it, the day to day variattion in human activity is far greater than the covid vs. post-covid differrences. Although we cannot be sure because about pre-covid human mobility, the available Google Mobility data indicate that human mobility within the parks might have not changed much during COVID and hence that the relevance of our study might be compromised.
## Google Mobility vs stringency
ggplot(s, aes(x = StringencyIndex, y = parks_percent_change_from_baseline, col = Country)) +
stat_smooth(method = 'lm', se = FALSE)+
stat_cor(method="pearson", size =2)+
geom_point()
We are not sure whether the test between FID and Google Mobility is meaningful because it looks at whether FID is a truly plastic trait that changes according to daily changes in human mobility. Nevertheless, I have tested for that, in general and for each country separately. Also, as raw data indicated that there might be a quadratic effect, I spedified also quadratic models and models that use only negative Google Mobility index and only positive Google Mobility index values. ### quick and dirty exploration with ggplot
ss[, Nsp := .N, by ='sp']
ss[, weekday:=factor(weekday, levels = c('Monday','Tuesday','Wednesday','Thursday','Friday','Saturday','Sunday'))]
ggplot(ss, aes(x = parks_percent_change_from_baseline, y = FID, col = Country, lty = year_, pch = year_)) +
stat_smooth(method='rlm', lwd =0.5) +
facet_wrap(~weekday, nrow =1)
ggplot(ss[Nsp>5], aes(y = FID, x = parks_percent_change_from_baseline, groups = sp, col = Country)) +
stat_smooth(method = 'lm', se = FALSE) +
labs(subtitle = "regresions for specis with >5 data points")
#theme(legend.position = 'none')
# geom_point(size =0.5, pch = 1) +
#facet_wrap(~sp, scales ='free_y')
# PREDICTIONS
# full model
ss[, country_weekday:=paste(Country,weekday)]
ms=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(scale(parks_percent_change_from_baseline)|genus)+ (1|Species)+(1|weekday) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc),
data = ss, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_ms = est_out(ms, 'ALL: (scale(google_mob)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (scale(google_mob)|Country) + (1|IDLocality) +(1|sp_loc)')
est_ms[, control_for_starting_distance:='yes']
mx=lmer(scale(log(FID))~
scale(Year)+
#scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(scale(parks_percent_change_from_baseline)|genus)+ (1|Species)+(1|weekday) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc),
data = ss, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_mx = est_out(mx, 'ALL: (scale(google_mob)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (scale(google_mob)|Country) + (1|IDLocality) +(1|sp_loc)')
est_mx[, control_for_starting_distance:='no']
msp=lmer(scale(log(FID))~ # singular due to weekday, else all ok
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(poly(parks_percent_change_from_baseline,2)|genus)+ (1|Species)+(1|weekday)+(1|sp_day_year) + (poly(parks_percent_change_from_baseline,2)|Country) + (1|IDLocality) +(1|sp_loc),
data = ss, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_msp = est_out(msp, 'ALL: (poly(google_mob,2)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (poly(google_mob,2)|Country) + (1|IDLocality) +(1|sp_loc)')
est_msp[, control_for_starting_distance:='yes']
mxp=lmer(scale(log(FID))~ # singular due to weekday, else all ok
scale(Year)+
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(poly(parks_percent_change_from_baseline,2)|genus)+ (1|Species)+(1|weekday)+(1|sp_day_year) + (poly(parks_percent_change_from_baseline,2)|Country) + (1|IDLocality) +(1|sp_loc),
data = ss, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_mxp = est_out(mxp, 'ALL: (poly(google_mob,2)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (poly(google_mob,2)|Country) + (1|IDLocality) +(1|sp_loc)')
est_mxp[, control_for_starting_distance:='no']
# PL
ps=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday)+(1|Species)+(1|sp_day_year),
data = ss[Country=='Poland'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_ps = est_out(ps, 'Poland: (1|weekday)+(1|Species)+(1|sp_day_year)')
est_ps[, control_for_starting_distance:='yes']
px=lmer(scale(log(FID))~
scale(Year)+
#scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday)+(1|Species)+(1|sp_day_year),
data = ss[Country=='Poland'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_px = est_out(px, 'Poland: (1|weekday)+(1|Species)+(1|sp_day_year)')
est_px[, control_for_starting_distance:='no']
psp=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday)+(1|Species)+(1|sp_day_year),
data = ss[Country=='Poland'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_psp = est_out(psp, 'Poland - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)')
est_psp[, control_for_starting_distance:='yes']
pxp=lmer(scale(log(FID))~
scale(Year)+
#scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday)+(1|Species)+(1|sp_day_year),
data = ss[Country=='Poland'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_pxp = est_out(pxp, 'Poland - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)')
est_pxp[, control_for_starting_distance:='no']
# HU
hs=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Hungary'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_hs = est_out(hs, 'Hungary: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_hs[, control_for_starting_distance:='yes']
hx=lmer(scale(log(FID))~
scale(Year)+
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Hungary'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_hx = est_out(hx, 'Hungary: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_hx[, control_for_starting_distance:='no']
hsp=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Hungary'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_hsp = est_out(hsp, 'Hungary - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_hsp[, control_for_starting_distance:='yes']
hxp=lmer(scale(log(FID))~
scale(Year)+
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Hungary'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_hxp = est_out(hxp, 'Hungary - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_hxp[, control_for_starting_distance:='no']
# CZ
cs=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Czech Republic'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_cs = est_out(cs, 'Czechia: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_cs[, control_for_starting_distance:='yes']
cx=lmer(scale(log(FID))~
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Czech Republic'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_cx = est_out(cx, 'Czechia: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_cx[, control_for_starting_distance:='no']
csp=lmer(scale(log(FID))~
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Czech Republic'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_csp = est_out(csp, 'Czechia - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_csp[, control_for_starting_distance:='yes']
cxp=lmer(scale(log(FID))~
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Czech Republic'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_cxp = est_out(cxp, 'Czechia - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_cxp[, control_for_starting_distance:='no']
# FI
fs=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Finland'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_fs = est_out(fs, 'Finland: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_fs[, control_for_starting_distance:='yes']
fx=lmer(scale(log(FID))~
scale(Year)+
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Finland'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_fx = est_out(fx, 'Finland: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_fx[, control_for_starting_distance:='no']
fsp=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Finland'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_fsp = est_out(fsp, 'Finland - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_fsp[, control_for_starting_distance:='yes']
fxp=lmer(scale(log(FID))~
scale(Year)+
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Finland'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_fxp = est_out(fxp, 'Finland - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_fxp[, control_for_starting_distance:='no']
# AU
as=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Australia'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_as = est_out(as, 'Australia: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_as[, control_for_starting_distance:='yes']
ax=lmer(scale(log(FID))~
scale(Year)+
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Australia'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_ax = est_out(ax, 'Australia: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_ax[, control_for_starting_distance:='no']
asp=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Australia'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_asp = est_out(asp, 'Australia - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_asp[, control_for_starting_distance:='yes']
axp=lmer(scale(log(FID))~
scale(Year)+
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|weekday) + (1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ss[Country=='Australia'], REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_axp = est_out(axp, 'Australia - quadratic: (1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality)')
est_axp[, control_for_starting_distance:='no']
# combine
o = rbind(est_ms, est_mx, est_msp, est_mxp,
est_as, est_ax, est_asp, est_axp,
est_cs, est_cx, est_csp, est_cxp,
est_hs, est_hx, est_hsp, est_hxp,
est_ps, est_px, est_psp, est_pxp,
est_fs, est_fx, est_fsp, est_fxp)
save(o, file = here::here('Data/dat_est_rev.Rdata'))
AIC(hs,hsp)
## df AIC
## hs 14 1959.335
## hsp 15 1954.779
AIC(cs,csp)
## df AIC
## cs 13 2500.136
## csp 14 2499.615
AIC(ps,psp)
## df AIC
## ps 13 742.5854
## psp 14 743.0644
AIC(as,asp)
## df AIC
## as 14 2456.332
## asp 15 2456.783
AIC(fs,fsp)
## df AIC
## fs 14 858.1777
## fsp 15 856.9257
Quadratic (indicated with ‘p’ in the model name) is never better than linear. ### PLOT estimates
load(here::here('Data/dat_est_rev.Rdata'))
o[predictor%in%c('scale(poly(parks_percent_change_from_baseline, 2))1','scale(parks_percent_change_from_baseline)'), predictor:='google mobility\n(linear)']
o[predictor%in%'scale(poly(parks_percent_change_from_baseline, 2))2', predictor:='google mobility\n(quadratic)']
oo = o[predictor %in% c('google mobility\n(linear)', 'google mobility\n(quadratic)')]
oo[, predictor := factor(predictor, levels = rev(c('google mobility\n(linear)', 'google mobility\n(quadratic)')))]
ggplot(oo, aes(x = estimate, y = predictor, color = model, shape = control_for_starting_distance)) +
geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
#geom_point(position = ggstance::position_dodgev(.6)) +
geom_point(position = position_dodge(width = width_), bg = 'white', size = 1.1) +
#scale_color_viridis(discrete=TRUE, begin=0, end = 0.5) +
#scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
#geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
# annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+
scale_shape_manual(name = 'Controlled for starting distance', guide = guide_legend(reverse = TRUE), values = c(23,21))+
scale_color_manual(name = "Model", guide = guide_legend(reverse = TRUE), values = col_)+
scale_x_continuous(breaks = round(seq(-0.6, 0.9, by = 0.1),1))+
ylab("") +
xlab("Standardized effect size\n[on flight initiation distance]") +
#coord_cartesian(xlim = c(-.15, .15)) +
#scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
theme_bw() +
theme(
legend.position = "right",
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
# legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
legend.margin=margin(0,0,0,0),
#legend.position=c(0.5,1.6),
plot.title = element_text(color = "grey", size = 7),
plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
# axis.text.x = element_text()
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(, size = 6),
axis.text.y = element_text(colour = "black", size = 7),
axis.title = element_text(size = 7)
)
load(here::here('Data/dat_est_rev.Rdata'))
o[predictor%in%c('poly(parks_percent_change_from_baseline, 2)1','scale(parks_percent_change_from_baseline)'), predictor:='google mobility\n(linear)']
oo = o[predictor %in% c('google mobility\n(linear)') & !model%in%o[predictor%in%'poly(parks_percent_change_from_baseline, 2)2', unique(model)]]
ggplot(oo, aes(x = estimate, y = model, shape = control_for_starting_distance)) +
geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
#geom_point(position = ggstance::position_dodgev(.6)) +
geom_point(position = position_dodge(width = width_), bg = 'white', size = 1.1) +
#scale_color_viridis(discrete=TRUE, begin=0, end = 0.5) +
#scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
#geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
# annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+
scale_shape_manual(name = 'Controlled for starting distance', guide = guide_legend(reverse = TRUE), values = c(23,21))+
# scale_color_manual(name = "Model", guide = guide_legend(reverse = TRUE), values = col_)+
ylab("") +
xlab("Standardized effect size\n[on flight initiation distance]") +
#coord_cartesian(xlim = c(-.15, .15)) +
#scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
theme_bw() +
theme(
legend.position = "right",
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
# legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
legend.margin=margin(0,0,0,0),
#legend.position=c(0.5,1.6),
plot.title = element_text(color = "grey", size = 7),
plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
# axis.text.x = element_text()
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(, size = 6),
axis.text.y = element_text(colour = "black", size = 7),
axis.title = element_text(size = 7)
)
# generate predictions for each country
nsim_= 5000
ss[, sin_ :=sin(rad)]
ss[, cos_ :=cos(rad)]
ss[, poly_lin := poly(parks_percent_change_from_baseline,2)[,1]]
ss[, poly_qua := poly(parks_percent_change_from_baseline,2)[,2]]
# HU, FI
p1 = foreach(j = c( 'Finland','Hungary'), .combine = rbind) %do% {
# j = 'Finland'
ssp = ss[Country == j]
psp=lmer(FID_ln~
Year+
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua +
(1|weekday)+(1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ssp, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
bsim = sim(psp, n.sim=nsim_)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5)) # coefficients
newD=data.frame(
Year = mean(ssp$Year),
SD_ln = mean(ssp$SD_ln),
flock_ln = mean(ssp$flock_ln),
body_ln = mean(ssp$body_ln),
rad = mean(ssp$rad),
Temp = mean(ssp$Temp),
parks_percent_change_from_baseline = seq(min(ssp$parks_percent_change_from_baseline), max(ssp$parks_percent_change_from_baseline), length.out = 100)) # values to predict for
newD$sin_ = sin(newD$rad)
newD$cos_ = cos(newD$rad)
newD$poly_lin = poly(newD$parks_percent_change_from_baseline,2)[,1]
newD$poly_qua = poly(newD$parks_percent_change_from_baseline,2)[,2]
X <- model.matrix(~
Year+
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua
, data=newD) # exactly the model which was used has to be specified here
# calculate predicted values and creditability intervals
newD$pred <-X%*%v #newD$fit_b <- plogis(X%*%v) # in case on binomial scaleback
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(i in 1:nsim) predmatrix[,i] <- X%*%bsim@fixef[i,]
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
newD$Country = j
return(newD)
}
# Poland
ssp = ss[Country=='Poland']
psp=lmer(FID_ln~
Year+
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua +
(1|weekday)+(1|Species)+(1|sp_day_year),
data = ssp, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
bsim = sim(psp, n.sim=nsim_)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5)) # coefficients
newD=data.frame(
Year = mean(ssp$Year),
SD_ln = mean(ssp$SD_ln),
flock_ln = mean(ssp$flock_ln),
body_ln = mean(ssp$body_ln),
rad = mean(ssp$rad),
Temp = mean(ssp$Temp),
parks_percent_change_from_baseline = seq(min(ssp$parks_percent_change_from_baseline), max(ssp$parks_percent_change_from_baseline), length.out = 100)) # values to predict for
newD$sin_ = sin(newD$rad)
newD$cos_ = cos(newD$rad)
newD$poly_lin = poly(newD$parks_percent_change_from_baseline,2)[,1]
newD$poly_qua = poly(newD$parks_percent_change_from_baseline,2)[,2]
X <- model.matrix(~
Year+
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua
, data=newD) # exactly the model which was used has to be specified here
# calculate predicted values and creditability intervals
newD$pred <-X%*%v #newD$fit_b <- plogis(X%*%v) # in case on binomial scaleback
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(i in 1:nsim) predmatrix[,i] <- X%*%bsim@fixef[i,]
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
newD$Country = 'Poland'
p2 = newD
# CZ
ssp = ss[Country=='Czech Republic']
psp=lmer(FID_ln~
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua +
(1|weekday)+(1|Species)+(1|sp_day_year)+ (1|IDLocality),,
data = ssp, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
bsim = sim(psp, n.sim=nsim_)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5)) # coefficients
newD=data.frame(
SD_ln = mean(ssp$SD_ln),
flock_ln = mean(ssp$flock_ln),
body_ln = mean(ssp$body_ln),
rad = mean(ssp$rad),
Temp = mean(ssp$Temp),
parks_percent_change_from_baseline = seq(min(ssp$parks_percent_change_from_baseline), max(ssp$parks_percent_change_from_baseline), length.out = 100)) # values to predict for
newD$sin_ = sin(newD$rad)
newD$cos_ = cos(newD$rad)
newD$poly_lin = poly(newD$parks_percent_change_from_baseline,2)[,1]
newD$poly_qua = poly(newD$parks_percent_change_from_baseline,2)[,2]
X <- model.matrix(~
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua
, data=newD) # exactly the model which was used has to be specified here
# calculate predicted values and creditability intervals
newD$pred <-X%*%v #newD$fit_b <- plogis(X%*%v) # in case on binomial scaleback
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(i in 1:nsim) predmatrix[,i] <- X%*%bsim@fixef[i,]
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
newD$Country = 'Czech Republic'
newD$Year = 2021
p3= newD
# AU
ssp = ss[Country == 'Australia']
psp=lmer(FID_ln~
Year+
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua +
(1|Species)+(1|sp_day_year) + (1|IDLocality),
data = ssp, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
bsim = sim(psp, n.sim=nsim_)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5)) # coefficients
newD=data.frame(
Year = mean(ssp$Year),
SD_ln = mean(ssp$SD_ln),
flock_ln = mean(ssp$flock_ln),
body_ln = mean(ssp$body_ln),
rad = mean(ssp$rad),
Temp = mean(ssp$Temp),
parks_percent_change_from_baseline = seq(min(ssp$parks_percent_change_from_baseline), max(ssp$parks_percent_change_from_baseline), length.out = 100)) # values to predict for
newD$sin_ = sin(newD$rad)
newD$cos_ = cos(newD$rad)
newD$poly_lin = poly(newD$parks_percent_change_from_baseline,2)[,1]
newD$poly_qua = poly(newD$parks_percent_change_from_baseline,2)[,2]
X <- model.matrix(~
Year+
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua
, data=newD) # exactly the model which was used has to be specified here
# calculate predicted values and creditability intervals
newD$pred <-X%*%v #newD$fit_b <- plogis(X%*%v) # in case on binomial scaleback
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(i in 1:nsim) predmatrix[,i] <- X%*%bsim@fixef[i,]
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
newD$Country = 'Australia'
p4 = newD
# combine
p = data.table(rbind(p1,p2,p3,p4))
p [, FID:=exp(pred)]
p [, FID_lwr:=exp(lwr)]
p [, FID_upr:=exp(upr)]
p[, Google_mobility:=parks_percent_change_from_baseline]
# plot
ggplot(p, aes(x = Google_mobility, y = pred, col = Country, fill = Country)) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=Google_mobility, fill = Country), alpha = 0.2, show.legend = NA, colour = NA) +
geom_line(aes(Google_mobility, y = pred, col = Country)) +
#scale_fill_viridis_c(option = "plasma", name = "°C\nwhen\npredated") +
#scale_x_continuous(expand = c(0, 0), lim = c(30,60), name = "Mid-day temperature [°C]") +
scale_y_continuous( name = "ln(Flight initiation distance)") +
#coord_cartesian(xlim = c(30,61), clip = 'off') +
#labs(tag = "(c)") +
theme_MB +
theme(
legend.text=element_text(size=5),
legend.title=element_text(size=6, hjust = 0.5)
)
# ' ### Global quadratic model
# generate predictions foo full model
nsim_= 5000
ss[, sin_ :=sin(rad)]
ss[, cos_ :=cos(rad)]
ss[, poly_lin := poly(parks_percent_change_from_baseline,2)[,1]]
ss[, poly_qua := poly(parks_percent_change_from_baseline,2)[,2]]
ssp =ss
psp=lmer(FID_ln~
Year+
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua +
(poly(parks_percent_change_from_baseline,2)|genus)+ (1|Species)+(1|weekday)+(1|sp_day_year) + (poly(parks_percent_change_from_baseline,2)|Country) + (1|IDLocality) +(1|sp_loc),
data = ssp, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
bsim = sim(psp, n.sim=nsim_)
v = apply(bsim@fixef, 2, quantile, prob=c(0.5)) # coefficients
newD=data.frame(
Year = mean(ssp$Year),
SD_ln = mean(ssp$SD_ln),
flock_ln = mean(ssp$flock_ln),
body_ln = mean(ssp$body_ln),
rad = mean(ssp$rad),
Temp = mean(ssp$Temp),
parks_percent_change_from_baseline = seq(min(ssp$parks_percent_change_from_baseline), max(ssp$parks_percent_change_from_baseline), length.out = 100)) # values to predict for
newD$sin_ = sin(newD$rad)
newD$cos_ = cos(newD$rad)
newD$poly_lin = poly(newD$parks_percent_change_from_baseline,2)[,1]
newD$poly_qua = poly(newD$parks_percent_change_from_baseline,2)[,2]
X <- model.matrix(~
Year+
SD_ln +
flock_ln +
body_ln +
sin_ + cos_ +
Temp+
poly_lin + poly_qua
, data=newD) # exactly the model which was used has to be specified here
# calculate predicted values and creditability intervals
newD$pred <-X%*%v #newD$fit_b <- plogis(X%*%v) # in case on binomial scaleback
predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
for(i in 1:nsim) predmatrix[,i] <- X%*%bsim@fixef[i,]
newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
p= data.table(newD)
p [, FID:=exp(pred)]
p [, FID_lwr:=exp(lwr)]
p [, FID_upr:=exp(upr)]
p[, Google_mobility:=parks_percent_change_from_baseline]
# plot
ggplot(p, aes(x = Google_mobility, y = pred)) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=Google_mobility), alpha = 0.2, show.legend = NA, colour = NA) +
geom_line(aes(Google_mobility, y = pred)) +
#scale_fill_viridis_c(option = "plasma", name = "°C\nwhen\npredated") +
#scale_x_continuous(expand = c(0, 0), lim = c(30,60), name = "Mid-day temperature [°C]") +
scale_y_continuous( name = "ln(Flight initiation distance)") +
#coord_cartesian(xlim = c(30,61), clip = 'off') +
#labs(tag = "(c)") +
theme_MB
ggplot(p, aes(x = Google_mobility, y = FID)) +
geom_ribbon(aes(ymin=FID_lwr, ymax=FID_upr, x=Google_mobility), alpha = 0.2, show.legend = NA, colour = NA) +
geom_line(aes(Google_mobility, y = FID)) +
#scale_fill_viridis_c(option = "plasma", name = "°C\nwhen\npredated") +
#scale_x_continuous(expand = c(0, 0), lim = c(30,60), name = "Mid-day temperature [°C]") +
scale_y_continuous( name = "Flight initiation distance [m]") +
#coord_cartesian(xlim = c(30,61), clip = 'off') +
#labs(tag = "(c)") +
theme_MB
# predictions
ssc <- ss[, sp_country_N := .N, by = sp_country]
ssc = ssc[sp_country_N>=15]
ssc[, log_FID :=log(FID)]
ssc[, log_SD :=log(SD)]
nsim = 1000
vv = list()
u = foreach(i = unique(ssc$sp_country), .combine = rbind) %do% {
# i = "Luscinia_megarhynchos\nHungary" "Gallinula_tenebrosa\nAustralia" "Columba_livia\nPoland" #ssc$sp_country[1]
ssci = ssc[sp_country == i]
if(ssci$Country[1] =='Poland' | ssci$sp_country[1]%in%c( "Gallinula_tenebrosa\nAustralia","Luscinia_megarhynchos\nHungary")){
mi <- lmer(log_FID ~
log_SD +
parks_percent_change_from_baseline +
(1 | weekday), data = ssci, REML = FALSE)
}else{
mi = lmer(log_FID ~
log_SD+
parks_percent_change_from_baseline +
(1 | weekday) + (1 | IDLocality), data = ssci, REML = FALSE)
}
bsim = sim(mi, n.sim = nsim)
v = apply(bsim@fixef, 2, quantile, prob = c(0.5)) # coefficients
newD = data.frame(
log_SD = mean(ssci$log_SD),
parks_percent_change_from_baseline = seq(min(ssci$parks_percent_change_from_baseline), max(ssci$parks_percent_change_from_baseline), length.out = 100)
) # values to predict for
X <- model.matrix(~log_SD +parks_percent_change_from_baseline, data = newD) # exactly the model which was used has to be specified here
# calculate predicted values and creditability intervals
newD$pred <- (X %*% v)
predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
for (j in 1:nsim) {
predmatrix[, j] <- (X %*% bsim@fixef[j, ])
}
#predmatrix[predmatrix < 0] <- 0
newD$pred <- as.numeric(newD$pred[, 1])
newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
newD$sp_country = as.character(i)
newD$Species= as.character(ssci$sp[1])
newD$Country = as.character(ssci$Country[1])
#print(i)
v = data.frame(v)
names(v) = 'estimate'
v$predictor =c('intercept','log_sd','google_mobiliity')
v = data.table(v)
v$sp_country = as.character(i)
v$Species= as.character(ssci$sp[1])
v$Country = as.character(ssci$Country[1])
vv[[i]] = data.table(v)
return(newD)
}
u = data.table(u)
vvv = do.call(rbind, vv)
u[, FID:=exp(pred)]
ggplot(u, aes(x = parks_percent_change_from_baseline, y = FID, col = Country, group = sp_country)) +
geom_line() +
labs(x = 'Google mobility index\n[% change from baseline]', y = 'Escape distance [m]', title = 'For species with N ≥ 15') +
scale_x_continuous(breaks = seq(-80, 140, by =20)) +
theme_bw() +
theme(
legend.position = "right",
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
# legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
legend.margin = margin(0, 0, 0, 1),
# legend.position=c(0.5,1.6),
plot.title = element_text(color ="grey40", size = 8),
#plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
#axis.line.y = element_blank(),
#axis.ticks.y = element_blank(),
axis.ticks = element_line(colour = ax_lines, size = 0.25),
# axis.text.x = element_text()
axis.ticks.length = unit(1, "pt"),
# axis.text = element_text(, size = 6),
axis.text = element_text(colour = "black", size = 7),
axis.title = element_text(size = 8)
)
# predictions
ss[parks_percent_change_from_baseline<0, google := 'before_zero']
ss[parks_percent_change_from_baseline>0, google := 'after_zero']
ss[, sp_country_google:= paste(sp_country, google)]
ssg <- ss[, sp_country_google_N := .N, by = sp_country_google]
ssg = ssg[sp_country_google_N >= 15]
ssg[, log_FID := log(FID)]
ssg[, log_SD := log(SD)]
nsim = 1000
ba = list()
uu = foreach(i = unique(ssg$sp_country_google), .combine = rbind) %do% {
# i = "Luscinia_megarhynchos\nHungary" "Gallinula_tenebrosa\nAustralia" "Columba_livia\nPoland" #ssc$sp_country[1]
ssgi = ssg[sp_country_google == i]
if(ssgi$Country[1] =='Poland' | ssgi$sp_country[1]%in%c( "Gallinula_tenebrosa\nAustralia","Luscinia_megarhynchos\nHungary")){
mi <- lmer(log_FID ~
log_SD +
parks_percent_change_from_baseline +
(1 | weekday), data = ssgi, REML = FALSE)
}else{
mi = lmer(log_FID ~
log_SD+
parks_percent_change_from_baseline +
(1 | weekday) + (1 | IDLocality), data = ssgi, REML = FALSE)
}
bsim = sim(mi, n.sim = nsim)
v = apply(bsim@fixef, 2, quantile, prob = c(0.5)) # coefficients
newD = data.frame(
log_SD = mean(ssgi$log_SD),
parks_percent_change_from_baseline = seq(min(ssgi$parks_percent_change_from_baseline), max(ssgi$parks_percent_change_from_baseline), length.out = 100)
) # values to predict for
X <- model.matrix(~log_SD +parks_percent_change_from_baseline, data = newD) # exactly the model which was used has to be specified here
# calculate predicted values and creditability intervals
newD$pred <- (X %*% v)
predmatrix <- matrix(nrow = nrow(newD), ncol = nsim)
for (j in 1:nsim) {
predmatrix[, j] <- (X %*% bsim@fixef[j, ])
}
#predmatrix[predmatrix < 0] <- 0
newD$pred <- as.numeric(newD$pred[, 1])
newD$lwr <- apply(predmatrix, 1, quantile, prob = 0.025)
newD$upr <- apply(predmatrix, 1, quantile, prob = 0.975)
newD$sp_country_google = as.character(i)
newD$sp_country = as.character(ssgi$sp_country[1])
newD$Species = as.character(ssgi$sp[1])
newD$Country = as.character(ssgi$Country[1])
newD$google = as.character(ssgi$google[1])
#print(i)
v = data.frame(v)
names(v) = 'estimate'
v$predictor =c('intercept','log_sd','google_mobiliity')
v = data.table(v)
v$sp_country_google = as.character(i)
v$sp_country = as.character(i)
v$Species= as.character(ssci$sp[1])
v$Country = as.character(ssci$Country[1])
v$google = as.character(ssgi$google[1])
ba[[i]] = data.table(v)
#l[[i]] = newD
return(data.table(newD))
}
bav = do.call(rbind, ba)
uu[, FID:=exp(pred)]
# plotting
ggplot(uu, aes(x = parks_percent_change_from_baseline, y = FID, col = Country, group = sp_country_google)) +
geom_line() +
labs(x = 'Google mobility index\n[% change from baseline]', y = 'Escape distance [m]', title = 'For species with N ≥ 15 per period') +
scale_x_continuous(breaks = seq(-80, 140, by =20)) +
theme_bw() +
theme(
legend.position = "right",
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
# legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
legend.margin = margin(0, 0, 0, 1),
# legend.position=c(0.5,1.6),
plot.title = element_text(color ="grey40", size = 8),
#plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
#axis.line.y = element_blank(),
#axis.ticks.y = element_blank(),
axis.ticks = element_line(colour = ax_lines, size = 0.25),
# axis.text.x = element_text()
axis.ticks.length = unit(1, "pt"),
# axis.text = element_text(, size = 6),
axis.text = element_text(colour = "black", size = 7),
axis.title = element_text(size = 8)
)
# model before
ssb = ss[google == 'before_zero']
#summary(factor(ssb$Country))
mbs=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(scale(parks_percent_change_from_baseline)|genus)+ (1|Species)+(1|weekday) + (1|sp_day_year) + (1|Country) + (1|IDLocality) +(1|sp_loc),
data = ssb, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_mbs = est_out(mbs, 'Before zero: (scale(google_mob)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (1|Country) + (1|IDLocality) +(1|sp_loc)')
est_mbs[, control_for_starting_distance:='yes']
mbx=lmer(scale(log(FID))~
scale(Year)+
#scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(scale(parks_percent_change_from_baseline)|genus)+ (1|Species)+(1|weekday) + (1|sp_day_year) + (1|Country) + (1|IDLocality) +(1|sp_loc),
data = ssb, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_mbx = est_out(mbx, 'Before zero: (scale(google_mob)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (1|Country) + (1|IDLocality) +(1|sp_loc)')
est_mbx[, control_for_starting_distance:='no']
mbsp=lmer(scale(log(FID))~ # singular due to weekday, else all ok
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(scale(poly(parks_percent_change_from_baseline,2))|genus)+ (1|Species)+(1|weekday)+(1|sp_day_year) + (1|Country) + (1|IDLocality) +(1|sp_loc),
data = ssb, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_mbsp = est_out(mbsp, 'Before zero: (poly(google_mob,2)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (1|Country) + (1|IDLocality) +(1|sp_loc)')
est_mbsp[, control_for_starting_distance:='yes']
mbxp=lmer(scale(log(FID))~ # singular due to weekday, else all ok
scale(Year)+
#scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(scale(poly(parks_percent_change_from_baseline,2))|genus)+ (1|Species)+(1|weekday)+(1|sp_day_year) + (1|Country) + (1|IDLocality) +(1|sp_loc),
data = ssb, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_mbxp = est_out(mbxp, 'Before zero: (poly(google_mob,2)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (1|Country) + (1|IDLocality) +(1|sp_loc)')
est_mbxp[, control_for_starting_distance:='no']
# model after
ssa = ss[google == 'after_zero']
mas=lmer(scale(log(FID))~
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(scale(parks_percent_change_from_baseline)|genus)+ (1|Species)+(1|weekday) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc),
data = ssa, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_mas = est_out(mas, 'After zero: (scale(google_mob)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (scale(google_mob)|Country) + (1|IDLocality) +(1|sp_loc)')
est_mas[, control_for_starting_distance:='yes']
max=lmer(scale(log(FID))~
scale(Year)+
#scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(parks_percent_change_from_baseline)+
(scale(parks_percent_change_from_baseline)|genus)+ (1|Species)+(1|weekday) + (1|sp_day_year) + (scale(parks_percent_change_from_baseline)|Country) + (1|IDLocality) +(1|sp_loc),
data = ssa, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_max = est_out(max, 'After zero: (scale(google_mob)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (scale(google_mob)|Country) + (1|IDLocality) +(1|sp_loc)')
est_max[, control_for_starting_distance:='no']
masp=lmer(scale(log(FID))~ # singular due to weekday, else all ok
scale(Year)+
scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(scale(poly(parks_percent_change_from_baseline,2))|genus)+ (1|Species)+(1|weekday)+(1|sp_day_year) + (scale(poly(parks_percent_change_from_baseline,2))|Country) + (1|IDLocality) +(1|sp_loc),
data = ssa, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_masp = est_out(masp, 'After zero: (poly(google_mob,2)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (poly(google_mob,2)|Country) + (1|IDLocality) +(1|sp_loc)')
est_masp[, control_for_starting_distance:='yes']
maxp=lmer(scale(log(FID))~ # singular due to weekday, else all ok
scale(Year)+
# scale(log(SD))+
scale(log(FlockSize))+
scale(log(BodyMass))+
scale(sin(rad)) + scale(cos(rad)) +
#scale(Day)+
scale(Temp)+
scale(poly(parks_percent_change_from_baseline,2))+
(1|genus)+ (1|Species)+(1|weekday)+(1|sp_day_year) + (1|Country) + (1|IDLocality) +(1|sp_loc),
data = ssa, REML = FALSE,
control = lmerControl(
optimizer ='optimx', optCtrl=list(method='nlminb'))
)
est_maxp = est_out(maxp, 'After zero: (poly(google_mob,2)|genus)+(1|Species)+(1|weekday) +(1|sp_day_year) + (poly(google_mob,2)|Country) + (1|IDLocality) +(1|sp_loc)')
est_maxp[, control_for_starting_distance:='no']
# combine
oab = rbind(est_mbs, est_mbx, est_mbsp, est_mbxp,
est_mas, est_max, est_masp, est_maxp)
save(oab, file = here::here('Data/dat_est_bef-aft_rev.Rdata'))
# plot
load(here::here('Data/dat_est_bef-aft_rev.Rdata'))
oab[predictor%in%c('scale(poly(parks_percent_change_from_baseline, 2))1','scale(parks_percent_change_from_baseline)'), predictor:='google mobility\n(linear)']
oab[predictor%in%'scale(poly(parks_percent_change_from_baseline, 2))2', predictor:='google mobility\n(quadratic)']
ooab = oab[predictor %in% c('google mobility\n(linear)', 'google mobility\n(quadratic)')]
ooab[, predictor := factor(predictor, levels = rev(c('google mobility\n(linear)', 'google mobility\n(quadratic)')))]
ggplot(ooab, aes(x = estimate, y = predictor, color = model, shape = control_for_starting_distance)) +
geom_vline(xintercept = 0, color = "grey", linetype = "dotted") +
geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0, position = ggstance::position_dodgev(width_)) +
#geom_point(position = ggstance::position_dodgev(.6)) +
geom_point(position = position_dodge(width = width_), bg = 'white', size = 1.1) +
#scale_color_viridis(discrete=TRUE, begin=0, end = 0.5) +
#scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) +
#geom_text( aes(x = n_pos,label = N), vjust = 0, size = 1.75, position = ggstance::position_dodgev(width_))+ # 3 positions for 3 bars
# annotate("text", x=log10(3), y=85, label= "Used", col = "grey30", size = 2.5)+
scale_shape_manual(name = 'Controlled for starting distance', guide = guide_legend(reverse = TRUE), values = c(23,21))+
scale_color_manual(name = "Model", guide = guide_legend(reverse = TRUE), values = col_)+
scale_x_continuous(breaks = round(seq(-0.6, 0.9, by = 0.1),1))+
ylab("") +
xlab("Standardized effect size\n[on flight initiation distance]") +
#coord_cartesian(xlim = c(-.15, .15)) +
#scale_x_continuous(breaks = round(seq(-.15, .15, by = 0.05),2)) +
theme_bw() +
theme(
legend.position = "right",
legend.title = element_text(size = 7),
legend.text = element_text(size = 6),
# legend.spacing.y = unit(0.1, 'cm'),
legend.key.height = unit(0.5, "line"),
legend.margin=margin(0,0,0,0),
#legend.position=c(0.5,1.6),
plot.title = element_text(color = "grey", size = 7),
plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r = 0.5, unit = "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = ax_lines, size = 0.25),
axis.line.y = element_blank(),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = ax_lines, size = 0.25),
# axis.text.x = element_text()
axis.ticks.length = unit(1, "pt"),
axis.text.x = element_text(, size = 6),
axis.text.y = element_text(colour = "black", size = 7),
axis.title = element_text(size = 7)
)
ss[, NspC := .N, by ='sp_country']
ssc = ss[NspC>9]
ssc[, sp2 := gsub(" ", "\n", sp)]
ssc[, Google_mobility:=parks_percent_change_from_baseline]
# two rows labels
g2 =
ggplot(ssc, aes(x = Google_mobility, y = FID,)) +
#stat_smooth(method = 'rlm', se = FALSE, col = 'black', lwd = 0.5)+
geom_point(pch = 21, alpha = 0.7, aes(fill = Country), col = 'white') +
stat_smooth(se = FALSE, aes(col = Country), lwd = 0.5)+ # show_guide=TRUE
facet_wrap(~sp2, ncol = 6) +
scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = FALSE)) +
scale_color_viridis(discrete=TRUE,guide = guide_legend(reverse = FALSE)) +
scale_x_continuous("Google mobility index", expand = c(0, 0)) +
scale_y_continuous("Flight initiation distance [m]", expand = c(0, 0), trans = 'log10') +
#annotate("text", x = 1, y = 1, label = c(rep("", 52),"Observation"), hjust = -0.08, size = 1) +
#labs(title = "Species means per sampling location")+
#scale_colour_manual(values=c('grey60'))+
#scale_color_manual(name = 'try', values = c('LOESS smoothed = "grey60"'))+
theme_MB +
theme(
plot.title = element_text(size=7),
strip.background = element_blank(),
#strip.text.x = element_text(size = 4.5, color="grey30", margin=margin(1,1,1,1,"mm")),
#panel.spacing = unit(1, "mm"),
legend.position = c(1, 0.01),
legend.justification = c(1, 0),
legend.title = element_blank(),
#legend.spacing.y = unit(-0.78, "cm")
#legend.spacing.y = unit(0.02, "cm") use if LOESS smooth text as legend
legend.spacing.y = unit(-0.9, "cm")
)
gg2 <- ggplotGrob(g2) #gg$layout$name
ggx2 <- gtable_filter_remove(gg2, name = c('axis-b-2-9','axis-b-5-8'), trim = FALSE)# paste0("axis-b-", c(2, 4), "-9")
grid.draw(ggx2)